/*
 * Example_CalciteDissolution.h
 *
 *  Created on: 19 Oct 2011
 *      Author: allan
 */

// C++ includes
#include <iostream>
#include <iomanip>
using namespace std;

// GeoReact includes
#include <GeoReact.h>

void Example_CalciteDissolution()
{
	const double ssaCalcite  = 9.8 * 100.09 * 1.0E-4; // conversion from cm2/g to m2/mol
	
	GeoManager manager;
	
	manager.SpeciateAqueousSpecies({"H2O", "Calcite"});
	manager.SetMineralSpecies({"Calcite"});
	
	MineralKineticsInfo kCalcite("Calcite");
	kCalcite.AddMineralMechanism(1.54881662E-6, 23.5);                                      // neutral mechanism
	kCalcite.AddMineralMechanism(0.501187234, 14.4, {make_tuple("H[+]", "activity", 1.0)}); // acid mechanism
	kCalcite.SetSpecificSurfaceArea(ssaCalcite);
	
	manager.SetMineralKinetics("Calcite", kCalcite);
	
	GeoSystem system = manager.AssembleGeoSystem();
	
	cout << system << endl;
	
	GeoEquilibrator equilibrator(system);
	
	equilibrator.ImposeChargeBalance();
	equilibrator.ImposeSpeciesMass("H2O", 55.508);
	equilibrator.ImposeSpeciesMass("Ca[2+]", 1.0E-16);
	equilibrator.ImposeSpeciesMass("HCO3[-]", 1.0E-16);
	equilibrator.ImposeSpeciesMass("Calcite", 1.0);
	
	GeoState state(system);
	
	state.SetTemperature(70.0);
	state.SetPressure(1.0);
	
	equilibrator.Equilibrate(state);
	
	cout << state << endl;
	
	GeoReactor reactor(system);
	
	const double hour = 3600;
	const double day  = 24 * hour;
	
	double t = 0.0, tEnd = day;
	
	cout << setprecision(10) << fixed;

	while(true)
	{
		cout << setw(20) << left << t/hour;
		cout << setw(20) << left << state.Get("Moles", "Calcite") * 100.09;
		cout << setw(20) << left << state.Get("pH") << endl;
		
		if(t == tEnd) break;
		
		reactor.React(t, tEnd, state);
	}
}
